1. IntroductionSurface roughness elements such as ribs or cubes are widely adopted to enhance flow mixing and heat transfer in engineering applications, because they can significantly modify flow properties and turbulent structures. [1–5] The influence of surface roughness on flow is foundational and essential for further heat and mass transfer. Therefore, it is important to investigate the mechanism of turbulent flow over a rough surface.
However, it has been a controversial issue about the influence region of surface roughness on turbulent flow. Townsend's wall-similarity hypothesis [6] states that the effects of surface roughness on turbulent motions are only limited in the roughness sub-layer, which is about five times of the roughness height; in the outer layer the mean flows and turbulent fluctuations are similar for smooth and rough wall cases. This hypothesis is supported by a number of studies such as Orlandi et al., [7, 8] Flack et al., [9] Wu and Christensen, [10] and Castro. [11] On the contrary, there are also several experimental and numerical results contradicting the wall-similarity hypothesis. It has been found that the outer flow is also affected by surface roughness, such as Krogstad and Antonia [12] for rod roughness, Lee and Sung [13] for rod and cube roughness, and Volino et al. [14] for rod roughness. Moreover, the velocity gradient in the overlap region leads to a logarithmic velocity profile, expressed by [13]
 | (1) |
where
κ is the von Karman constant and
B is the additive intercept at
.
and
have been widely used for all the wall-bounded flows. However, some experimental results in the atmospheric surface layer [
15–
17] showed that
κ is not universal and decreases with the increase of roughness Reynolds number (
, where
is friction velocity and
y
0 is defined such that
at
. This was contradicted by the latter measurement in the atmospheric surface layer. [
18] Recently, Nagib and Chauhan [
19] suggested that
κ is different in boundary layer, pipes, and channels. Leonardi
et al. [
20,
21] found that
κ varies in ribbed- and cubic-channel flow by direct numerical simulation (DNS) study. It was also found that
κ is changed in open-channel flows with bed-load. [
22,
23] Therefore, the influences of surface roughness on turbulent flow have been insufficiently understood and need to be further investigated.
As for the research methods, it has been quite challenging to measure accurately near-wall in experiments, as the flow over rough surface is very complex. [24] Therefore, numerical methods such as DNS and large eddy simulation (LES) are employed in most of the research. Leonardi et al. [25] compared DNS and LES methods coupled with immersed boundary method (IBM) to study turbulent channel flow over ribbed-rough wall, and found that DNS and LES can both resolve such rough-wall turbulence very well.
To investigate the influences of surface cube roughness on turbulent velocity profiles as well as coherent structures, here we adopt LES together with IBM to simulate turbulent channel flow over a cube-rough wall. Meanwhile, the wall-similarity hypothesis is further examined for different roughness height.
2. Simulation methodThe computational channel with surface cubes is shown in Fig. 1. The channel size is
in the x,y, and z directions, respectively. The criterion of the minimal channel width is that it should be wide enough to contain an ejection and its associated sweep. According to [26] and [27], the size of
is sufficient for investigating the turbulent channel flow. For the rough wall cases, three-dimensional cube array is arranged on the bottom wall of the channel with equal streamwise and spanwise spacing, while the upper channel wall is smooth. The calculation cases in detail are shown in Table 1. The roughness height is
, or
, respectively. The length and width of the roughness are both
. The roughness spacings w are also
in the streamwise and spanwise directions for all the rough wall cases.
Table 1.
Table 1.
 Table 1.
Computational cases.
.
| Case No. |
Average roughness |
Roughness height |
Roughness |
|
height k/h
|
in wall units
|
spacing w/h
|
| A |
0 |
0 |
0 |
| B |
0.05 |
9 |
0.375 |
| C |
0.1 |
18 |
0.375 |
| D |
0.15 |
27 |
0.375 |
| Table 1.
Computational cases.
. |
The friction velocity
is determined by the wall shear stress
and the fluid density ρ, calculated by
 | (2) |
For rough-channel flow cases, the wall shear stress on the rough wall
is different from that on the smooth wall
. In this study, the flow is driven by a constant pressure drop
. Therefore, the driven force is balanced with the sum of smooth and rough wall shear stresses, that is
 | (3) |
On the top wall (smooth), the wall shear stress
can be calculated as
 | (4) |
With Eqs. (
4) and (
4) the wall shear stress of bottom rough wall
can be calculated as
 | (5) |
Then the friction velocity for the top wall
and that for the bottom wall
can be computed by the Eq. (
2) with their respective wall shear stress.
For rough-channel flow case, the friction Reynolds number is defined by the global mean friction velocity
, which is calculated from pressure gradient
. The friction Reynolds number is computed by
 | (6) |
which is 180 in the study. Therefore the roughness height
is
, and 27 in wall units.
LES coupled with IBM are used to simulate the channel flow with cube surface roughness. The governing equations after a classical box filter are written as
 | (7) |
 | (8) |
where
is the velocity component,
p is pressure,
f
i
is the momentum forcing defined at the roughness element faces inside the body.
is the constant pressure drop to drive the channel flow, which is 1 in the streamwise direction and 0 in the spanwise and wall-normal directions.
is the sub-grid-scale stress, which is approximated by the eddy viscosity hypothesis
 | (9) |
The eddy viscosity
and strain rate tensor
are defined as
 | (10) |
 | (11) |
where
is an empirical constant and taken as 0.1 here. For the channel flow, a damping function
is applied to Eq. (
10), because the Smagorinsky model has overlarge dissipation predictions near the wall and the viscosity must be zero at the wall. This damping function was also applied in [
28] for rough wall flows.
Equations (7) and (8)are solved using a fractional step method by Chorin et al. [29] The viscous term is discretized by second-order central scheme and the advection term is discretized using second-order hybrid scheme. The momentum equations are integrated explicitly using the third-order Runge–Kutta algorithm. The pressure Poisson equation is solved by means of Fourier series expansions in the streamwise and spanwise directions with tridiagonal matrix inversion.
The roughness elements are identified by IBM, which is first proposed by Peskin. [30] In the present study, the direct discrete forcing approach of IBM [31] is adopted to model the surface roughness. To ensure the non-slip boundary condition on the solid wall of the surface roughness, the direct forcing f
i
is imposed in the momentum equation, that is
 | (12) |
where
and
are the velocities for time step
n and the next time step
contains advection term, pressure gradient term, and diffusion term.
f
i
is obtained by
 | (13) |
where
is obtained by using interpolation method at the roughness element surfaces
.
f
i
is zero for the non-solid region. The detailed description of the direct discrete forcing approach of IBM can be found in [
30].
Periodic boundary conditions are applied in the streamwise and spanwise directions. At the initial time, the parabolic velocity distribution overlaid random velocity fluctuation are imposed in the inlet of the channel, and then the flow develops to turbulence state with the computational time. To calculate the rough wall cases, 1.04 million staggered grids are used in total, shown in Fig. 2(a). The grids are uniform in the streamwise and spanwise directions, while grid points are clustered near-wall in the wall-normal direction using a hyperbolic tangent function. Figure 2(b) displays the enlarged view of near-wall grids.
In the wall-normal direction, the roughness heights are 9, 18, and 27, respectively. The first grid spacing around the roughness elements
in wall units, and the maximum grid spacing on the crest of the roughness elements
. Thus there are
, and 20 grids in the wall-normal direction for the three rough wall cases, respectively. In the streamwise and spanwise directions, the grids are uniform and
in wall units. The length and width of the roughness are both
. Therefore, there are 8 grids in the streamwise and spanwise directions for roughness elements. The Kolmogorov scale of turbulent channel flow at the present Reynolds number is about 2 wall units. [32] Thus the present grid space never exceeds 5 times of the Kolmogorov scale around the roughness elements. Therefore, the grid resolution is fine enough for LES.
When the computation sustains several flow periods from the initial condition, the turbulent flow will become fully developed. In the practical simulation, we monitored the dimensionless instantaneous velocity at the channel center versus dimensionless computational time, as shown in Fig. 3. If the dimensionless instantaneous velocity curve is fully stable over time, it indicates that the calculated turbulent flow has been fully developed. The mean velocity and velocity fluctuation statistics are conducted after more than 100 channel flow periods of velocity convergence, when the turbulent flow velocity and the Reynolds stress already reach the steady state.
To verify the present codes of LES and IBM, the mean and fluctuating velocity profiles for turbulent channel flow over smooth and rough walls are compared with literature data, shown in Figs. 4 and 5. For smooth wall case, the results for both mean and fluctuating velocity profiles agree very well with DNS data by Kim et al. [32] For rough wall case, turbulent channel flow over ribbed surface is calculated according to the experimental case by Burattini et al. [33]
The grid independent check was conducted as part of validation, shown in Fig. 5. The grid numbers are
and
. The grid spacings of coarse mesh near the rough wall are
,
near the wall and
at the center of the channel. For the fine mesh, the grid spacings near the rough wall are
,
and
. From Fig. 5, it can be found that there are some differences both on the mean and fluctuating velocities between the numerical results of coarse mesh and experimental data, while the results of the fine mesh are in good agreement with the experimental data. This indicates that the present fine grids can resolve turbulent channel flow with rough walls very well. More details about the flow configuration of the ribbed-wall flow can be found in [33]. Therefore, the present codes are able to solve accurately turbulent channel flow both for smooth and rough walls.
3. Results and discussion3.1. Turbulent velocity profilesThe mean velocity profiles of turbulent flows over smooth and rough walls are shown in Fig. 6. For rough wall cases, the friction velocities are different on the top (smooth) wall from the value on the bottom (rough) wall. Thus in Figs. 6–9, the velocity and stress profiles are normalized by
for the smooth wall cases and by
for the rough wall cases. These normalizations are also consistent with [33].
It can be seen that the mean velocities for the rough wall cases are all decreased compared with the smooth wall case for the extra form drag induced by the roughness elements. Moreover, the mean velocity is decreased with the increase of the roughness height, because higher roughness elements will produce more form drag.
For the smooth-wall flow, a perfect logarithmic velocity distribution can be observed in Fig. 6, and can be expressed as Eq. (1), where
and
. For the rough-wall flows, there are still logarithmic regions for the mean velocity profiles. This indicates that the log-law of the wall is not destroyed by surface cube array with different roughness height. Further, the velocity distribution in the logarithmic region for the rough wall cases can be expressed as
 | (14) |
where
is called the roughness function representing the increase of flow drag induced by roughness. It is equal to 0 for smooth wall case A and is 0.6, 1.7, and 2.5 for the rough wall cases B, C, and D, respectively. Therefore, it can be concluded that the roughness function
increases obviously with the roughness height.
Moreover, it can be found that the slopes of the mean velocity profiles in the logarithmic region are significantly different between smooth wall and rough wall cases. The slope is increased with the increase of roughness height
, which means that the von Karman constant decreases with the increase of roughness height. That is consistent with [21]. As the surface cubes probably modify the turbulent coherent structures in the wall-bounded flow, the balance between the production and dissipation of turbulent kinetic energy is broken in the rough-wall channel flow. This may cause the variation of von Karman constant. [21] Therefore, the present LES results support that the von Karman constant is decreased in the rough-wall channel flows, at least in the cases of this study.
The defect profiles of mean velocities for rough-wall channel flows are shown in Fig. 7. It can be seen that the mean velocity profiles for rough-wall flows are defective in the inner layer compared with the smooth wall case. This tendency becomes more significant with the increase of roughness height
. In the outer layer, the differences on the mean velocity profiles between the rough wall and smooth wall cases are relatively small. It can be concluded that the mean velocity profile is similar in the outer layer for smooth and rough wall cases.
The fluctuating velocity profiles of turbulent flows over smooth and rough walls are shown in Fig. 8. The critical positions between roughness sub-layer and outer layer for the three rough-wall cases are marked out. It can be observed that the streamwise fluctuating velocity is obviously decreased while the spanwise and wall-normal fluctuating velocities are significantly increased by the surface cube array in roughness sub-layer (
). These influences are highly associated with the modification of turbulent structures by surface roughness near the wall, which will be discussed later. Nevertheless, the turbulent fluctuating velocities are almost not affected by surface roughness in the outer layer (
) for the rough wall cases B and C, compared with the smooth wall case. For the rough wall case D, it seems that the turbulent fluctuations are influenced by surface roughness in even half of the channel shown in Figs. 8(a)–8(c). However, noting that the roughness height for case D is as high as
, the roughness sub-layer of case D is thus very thick to the center of the channel (
). It can be seen that the values of fluctuating velocities for case D return yet to approximate those for smooth wall case at y = h. This indicates that the wall-similarity hypothesis still holds even for high roughness. In consequence, these results support the wall-similarity hypothesis of Townsend for cube-rough wall flows.
The distributions of Reynolds shear stresses
for smooth and rough walls are shown in Fig. 9. In the roughness sub-layer (
), the Reynolds shear stress is significantly increased with the increase of roughness height. As Reynolds shear stress represents the mean momentum flux due to the fluctuating velocity, larger shear stress indicates that more intense momentum transport and the exchange are induced by surface roughness. However, the Reynolds shear stresses in the outer layer (
) are also not affected dramatically by surface roughness for all the rough wall cases, which is similar to the case of the velocity fluctuations. It proves again that the Townsend's wall-similarity hypothesis holds for the present surface roughness.
3.2. Coherent structuresQuasi-streamwise vortices (QSV) are typical coherent structures near the wall, which are highly associated with turbulent production and energy transport in boundary layer. To investigate the influence of the surface roughness array on QSV, instantaneous streamwise vorticity
is visualized numerically in the X–Y plane at
for smooth and rough walls, shown in Fig. 10. Moreover, time- and streamwise-averaged QSV (
) for smooth and rough wall cases are obtained to investigate the effects of surface roughness on the spatio-temporal averaged QSV, as shown in Fig. 11. The streamwise vorticity
is defined as
 | (15) |
In Figs. 10 and 11, the red and blue patterns in the near-wall regions represent the QSV pairs. The colors in the figures represent the magnitude of the vorticity. Large-scale QSV can also be observed for smooth and rough wall cases. This indicates that near-wall instantaneous QSV still exist above the roughness array. However, there is significant difference between smooth and rough wall cases. The generation positions, the sizes and the quantities of QSV are all modified by the roughness array. Compared with the smooth wall case, more intense and larger number of QSV pairs can be found for rough wall cases. Moreover, more intense ejection motions are induced from the cavities of the roughness array by larger size of QSV. The strength of ejection motion above the roughness array is also increased with the increase of roughness height. This is because that the turbulent kinetic energy is transported from the streamwise direction to the spanwise and wall-normal directions. This will enhance turbulent velocity fluctuations and momentum transport in wall-normal and spanwise directions. It is in agreement with the results in Fig. 7. Therefore, the turbulence is modified to be more isotropic by the introduction of the roughness array.
